##
## Replication file for main article analysis of Brad Epperly, Chris Witko, Ryan Strickler, and Paul White, "Rule by Violence, Rule by Law," in Perspectives on Politics
## See readme file and appendix replication file for additional information
##

setwd("~/fillintheblank")

library(ggplot2) # Figures
library(extrafont) # Fonts
library(lme4) # Multilevel models
library(simcf) # For creating subsets of data, see Chris Adolph's webpage if it's not available on CRAN
library(memisc) # Tables
library(texreg) # Tables
options(scipen=120,digits=2) # Output format

# Below, change to wherever you have downloaded the data:
alldata2 <- read.csv("~/Desktop/epperly_etal_lynching_perspectives_data.csv")

#
# Table 1 is just dates of policies, no replication via R
#

#
# Some code to reproduce Figure 2, which plots lynchings over time:
#
overtime <- tapply(alldata2$confirmed_black_lynching,alldata2$year,sum)			
length(overtime)
summary(alldata2$year)
years <- seq(1876,1952,1)
length(years)

plot(years,overtime ,type="l") # Base R plot, just to get an impression

tapply(alldata2$year[alldata2$fulladditivejimcrow > 0],alldata2$State[alldata2$fulladditivejimcrow > 0],min,na.rm=TRUE) # The year (after) the first lynching policy is implemented in each state
tapply(alldata2$year[alldata2$fulladditivejimcrow > 1],alldata2$State[alldata2$fulladditivejimcrow > 1],min,na.rm=TRUE) # The year (after) the second lynchinig policy is implemented in each state

# Create dataframe of relevant dates and yearly counts for ggplotting
dat1 <- cbind(years,overtime)
dat1 <- as.data.frame(dat1)

plot <- ggplot(dat1, aes(x = years))
plot <- plot + geom_line(aes(y = overtime), linetype = 1, size = 2)
plot <- plot + geom_vline(xintercept=mean(c(1893,1891,1889,1877,1891,1897,1890,1899,1882,1889,1894)),linetype=2,size=1.25)
plot <- plot + geom_vline(xintercept=mean(c( 1893, 1892, 1889, 1900, 1897, 1888, 1899, 1894, 1890, 1902  )),linetype=2,size=1.25)
plot <- plot + scale_x_continuous(name='Year')
plot <- plot + scale_y_continuous('Number of lynching events')
plot <- plot + theme_bw()
plot <- plot + theme(plot.title=element_text(size=25))
plot <- plot + theme(axis.title.y = element_text(size=16,angle=90,family="Helvetica", face="bold"),
          axis.text.y=element_text(size=16,family="Helvetica", face="bold"))
plot <- plot + theme(axis.title.x = element_text(size=16,vjust=0,family="Helvetica", face="bold"),
          axis.text.x  = element_text(size=16,family="Helvetica", face="bold"))
plot
#ggsave(filename="~/Dropbox/lynchingnew/POPsub/figure1.pdf",width=8,height=4,units="in")		

## 
## How to make Table 2:
## 
form1 <- confirmed_black_lynching ~  daystoelection1 + pr + scale(ppop) + 
			scale(percentblack100) + scale(I(percentblack100^2))  + 
			 scale(cotratio) + scale(year) + State + num_confirmed_black_lynchings
# Want to create the datasets for the pre/post Jim Crow data:
pre_anytwo <- subset(alldata2, fulladditivejimcrow <2)
post_anytwo <- subset(alldata2, fulladditivejimcrow>=2)
pre_anytwo <- extractdata(form1, pre_anytwo,na.rm=TRUE) # Just have the variables of interest, remove missing data
post_anytwo <- extractdata(form1, post_anytwo,na.rm=TRUE) # Just have the variables of interest, remove missing data

# Models for Table 2, with covariates centered and scaled
fe.twopre1 <- glm(confirmed_black_lynching ~  scale(daystoelection1) + scale(ppop) + scale(pr) +  
			scale(percentblack100) + scale(I(percentblack100^2))  + 
			 scale(cotratio) + scale(year) + State,
			family=binomial(link="logit"),data= pre_anytwo)			
hlm.twopre1 <- glmer(confirmed_black_lynching ~  scale(daystoelection1) +scale(ppop) + scale(pr) + 
			scale(percentblack100) + scale(I(percentblack100^2))  + 
			 scale(cotratio) + scale(year) + (1|State),
			family=binomial(link="logit"),data= pre_anytwo)			
fe.twopost1 <- glm(confirmed_black_lynching ~  scale(daystoelection1) +scale(ppop) + scale(pr) + 
			scale(percentblack100) + scale(I(percentblack100^2))  + 
			 scale(cotratio) + scale(year) + State,
			family=binomial(link="logit"),data= post_anytwo)			
hlm.twopost1 <- glmer(confirmed_black_lynching ~  scale(daystoelection1) +scale(ppop) + scale(pr) + 
			scale(percentblack100) + scale(I(percentblack100^2))  + 
			 scale(cotratio) + scale(year) + (1|State),
			family=binomial(link="logit"),data= post_anytwo)			

screenreg(list(fe.twopre1, hlm.twopre1,fe.twopost1, hlm.twopost1),
	stars=c(0.001,0.01,0.05),digits=2,booktabs=TRUE,dcolumn=TRUE,
	custom.coef.names=c(
	"(Intercept)",
	"Days to election",
	"Percent Populist vote",
	"Percent Republican vote",
	"Percent Black",
	"Percent Black squared",	
	"Cotton dependence",
	"Year",
	"Arkansas",
	"Florida",
	"Georgia",
	"Kentucky",
	"Louisiana",
	"Mississippi",
	"North Carolina",
	"South Carolina",
	"Tennessee",
	"Virginia"
	))
# The state effects are commented out, so as to not be included in the manuscript, but included here for the curious; change screenreg to texreg for Latex output


#
# Making plots (a) and (b) in Figure 2, showing the marginal effects of the two political covariates pre- and post-Jim Crow
#
library(effects)

#plot (a)
# removing the standardization for days until election to more easily plot its effects because we do not have to mess with the x-axis -- substantive results are unchanged
z1 <- glm(confirmed_black_lynching ~  daystoelection1 +scale(ppop) + scale(pr) + 
			scale(percentblack100) + scale(I(percentblack100^2))  + 
			 scale(cotratio) + scale(year) + State,
			family=binomial(link="logit"),data= pre_anytwo)
# Can use the effects library to get marginal effects, but take the output and plot it in ggplot for better display
plot(effect("daystoelection1",z1,xlevels=706,
	confidence.level=0.95,se=TRUE))
e1 <- effect("daystoelection1",z1,xlevels=706)
dat <- data.frame(x = seq( min(pre_anytwo$daystoelection1), max(pre_anytwo$daystoelection1),length.out=706), lower = summary(e1)$lower, 
	middle = summary(e1)$effect, upper = summary(e1)$upper)
	
z2 <- glm(confirmed_black_lynching ~  daystoelection1 +scale(ppop) + scale(pr) + 
			scale(percentblack100) + scale(I(percentblack100^2))  + 
			 scale(cotratio) + scale(year) + State,
			family=binomial(link="logit"),data= post_anytwo)		
e2 <- effect("daystoelection1",z2,xlevels=706)
dat2 <- data.frame(x = seq( min(post_anytwo$daystoelection1), max(post_anytwo$daystoelection1),length.out=706), lower = summary(e2)$lower, 
	middle = summary(e2)$effect, upper = summary(e2)$upper)
	
plot <- ggplot(dat, aes(x = x))
plot <- plot + geom_line(aes(y = lower), linetype = 2)
plot <- plot + geom_line(aes(y = middle), linetype = 1,size=2)
plot <- plot + geom_line(aes(y = upper), linetype = 2)
plot <- plot + geom_line(data=dat2,aes(y = lower), linetype = 2,size=1,color="gray")
plot <- plot + geom_line(data=dat2,aes(y = middle), linetype = 1,size=2,color="gray")
plot <- plot + geom_line(data=dat2,aes(y = upper), linetype = 2,size=1,color="gray")
plot <- plot + scale_x_continuous(name='')
plot <- plot + scale_y_continuous('Pr(Lynching)',breaks=c(0.001,0.002,0.003,0.004))
plot <- plot + theme_bw()
plot <- plot + theme(plot.title=element_text(size=25))
plot <- plot + theme(axis.title.y = element_text(size=16,angle=90,		
	family="Helvetica",face="bold"),
          axis.text.y=element_text(size=16,family="Helvetica", face="bold"))
plot <- plot + theme(axis.title.x = element_text(size=16,vjust=0,family="Helvetica", face="bold"),
          axis.text.x  = element_text(size=16,family="Helvetica", face="bold"))
plot
#ggsave(filename="~/Dropbox/lynchingnew/POPsub/figure2_plot_a.pdf",width=4,height=4,units="in")	

# plot (b)
# removing the standardization for percent populist to more easily plot its effects because we do not have to mess with the x-axis -- substantive results are unchanged
z3 <- glm(confirmed_black_lynching ~  scale(daystoelection1) +scale(pr) + ppop + 
			scale(percentblack100) + scale(I(percentblack100^2))  + 
			 scale(cotratio) + scale(year) + State,
			family=binomial(link="logit"),data= pre_anytwo)	
z4 <- glm(confirmed_black_lynching ~  scale(daystoelection1) +scale(pr) + ppop + 
			scale(percentblack100) + scale(I(percentblack100^2))  + 
			 scale(cotratio) + scale(year) + State,
			family=binomial(link="logit"),data= post_anytwo)	
e3 <- effect("ppop",z3,xlevels=86)
e4 <- effect("ppop",z4,xlevels=79)
dat3 <- data.frame(x = seq(0,85,1), lower = summary(e3)$lower, 
	middle = summary(e3)$effect, upper = summary(e3)$upper)
dat4 <- data.frame(x = seq(0,78,1), lower = summary(e4)$lower, 
	middle = summary(e4)$effect, upper = summary(e4)$upper)
	
plot <- ggplot(dat3, aes(x = x))
plot <- plot + geom_line(aes(y = lower), linetype = 2,size=1)
plot <- plot + geom_line(aes(y = middle), linetype = 1,size=2)
plot <- plot + geom_line(aes(y = upper), linetype = 2,size=1)
plot <- plot + geom_line(data=dat4,aes(y = lower), linetype = 2,size=1,color="gray")
plot <- plot + geom_line(data=dat4,aes(y = middle), linetype = 1,size=2,color="gray")
plot <- plot + geom_line(data=dat4,aes(y = upper), linetype = 2,size=1,color="gray")
plot <- plot + scale_x_continuous(name='')
plot <- plot + scale_y_continuous('Pr(Lynching)',limits=c(0,0.013))
plot <- plot + theme_bw()
plot <- plot + theme(plot.title=element_text(size=25))
plot <- plot + theme(axis.title.y = element_text(size=16,angle=90,		
	family="Helvetica", face="bold"),
          axis.text.y=element_text(size=16,family="Helvetica", face="bold"))
plot <- plot + theme(axis.title.x = element_text(size=16,vjust=0,family="Helvetica", face="bold"),
          axis.text.x  = element_text(size=16,family="Helvetica", face="bold"))
plot
#ggsave(filename="~/Dropbox/lynchingnew/POPsub/figure2_plot_b.pdf",width=4,height=4,units="in")	
